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Abstract 

We study thermal behavior of a recently introduced Hartree ensemble approximation, 
which allows for non-perturbative inhomogeneous field configurations as well as for ap- 
proximate thermalization, in the y? 4 model in 1+1 dimensions. Using ensembles with a 
free field thermal distribution as out-of-equilibrium initial conditions we determine ther- 
malization time scales. The time scale for which the system stays in approximate quantum 
thermal equilibrium is an indication of the time scales for which the approximation method 
stays reasonable. This time scale turns out to be two orders of magnitude larger than the 
time scale for thermalization, in the range of couplings and temperatures studied. We 
also discuss simplifications of our method which are numerically more efficient and make 
a comparison with classical dynamics. 
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1 Introduction 



It is highly desirable to be able to numerically simulate quantum field dynamics in real 
time. This will give an important tool for the study of non-perturbative phenomena in out- 
of-equilibrium systems, such as phase transitions in the early universe or the quark-hadron 
transition in heavy-ion collisions. Using real time dynamics may also offer an alternative to 
simulating equilibrium physics, just like molecular dynamics simulations provide a fruitful 
alternative to Monte Carlo simulations in other areas of physics. Simulating quantum fields 
in real time is very difficult. Direct approaches, such as solving the Schrodinger equation 
for the field wave-functional or evaluating the Minkowski path integral using Monte Carlo 
methods, are prohibitively time-consuming. One has to resort to approximate methods, of 
which the classical approximation ( see e.g. [1-10]), the Hartree approximation and large 



n methods (see e.g. IO, I2|, [13]]) are most commonly used. 

In the classical approximation one assumes that the fields follow classical equations of 
motion. This is reasonable when the occupation numbers of the field quanta are large, 
but in field theory this is never the case for all modes. For instance, at high temperatures 
the low momentum modes of the fields are highly occupied and follow the classical Boltz- 
mann distribution, but at large momenta occupation numbers are low and the classical 
distribution differs significantly from the quantum Bose-Einstein distribution, giving rise 
to Rayleigh- Jeans divergences. 

These divergences are absent in the Hartree and large n approximations, which include 
quantum effects in the field dynamics. In these approximations the density operator is 
effectively gaussian, such that all information is contained in the mean field and two- 
point function. These are usually assumed to be translationally invariant. A problem 
which arises under these circumstances is that the system does not thermalize [|ll], [12| , 
in contrast to the classical approximation which has no such problem M. The particles 
corresponding to the quantum modes of the field interact via the mean field and since 
this field is homogeneous there is no scattering which leads to redistribution of occupation 
numbers over different momentum modes and there is no thermalization. One way to 
amend this situation may be an improved approximation which includes direct scattering 

g. 



Recently we have extended the Hartree approximation by writing an initial density 
operator as an ensemble of coherent states with generally inhomogeneous mean fields 
and two-point functions. Some further discussion may be required to understand the 
motivation for this approach. To start, we note that there is a class of density operators 



p which can be written as a superposition of gaussian pure states (see |L5], [16] and the 
appendix A of |T7[):Q 

p= [d<pdir]p[<p,ir]\<p,ir)(<p,ir\. (1) 



The \<p,7r) are coherent states centered around <p(x) = ((p,Tr\tp(x)\tp,ir) and tt(x) 
(tp, 7r|7r(x)|^, 7r), and p[<p, tc] is a functional representing the density operator p. 

§ Operators are indicated with a caret. 
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For example, for a free scalar field the canonical distribution p = exp(— /3H[<p, it]) is 
represented as (see the appendix A of [T7j for a derivation) 



p[<p, tt] oc H exp [-(e^ k - 1)(tt£ + ^)/2cu k ] , (2) 

k 

where k labels the modes of the field with frequency uy. 

By writing the density operator in this form, and we emphasize that e.g. for the canon- 
ical distribution (^) there are no approximations involved, we have achieved four things. 
Firstly, we have made contact with the classical approximation. If the mean field in 
a coherent state is large compared to the width of the state, the gaussian wavepacket 
approximately follows a classical trajectory and the mean field can be thought of as a 
classical field. This then suggests that the individual coherent states in the ensemble may 
be referred to as "realizations" . However, by using an ensemble of coherent states rather 
than classical fields, we have a (hopefully much) better description for those modes that 
have low occupation numbers for which the classical dynamics is a poor approximation. 
Secondly, we have expressed a (typically non-gaussian) initial density operator in terms 
of gaussian states. These are optimal for the Hartree method, which we want to use to 
approximate the dynamics of these states. Thirdly, the mean fields in the individual coher- 
ent states are inhomogeneous, therefore the particles can interact with the inhomogeneous 
mean field, such that the energy may get distributed over the full momentum range. In 



1 17] we found that this leads to approximate thermalization in coarse grained distributions. 
Averaging over the initial ensemble is not necessary per se, since it also occurs in each 
individual member of the ensemble, because of coarsening over the volume, provided the 
volume is large enough to contain a sufficient number of decorrelated systems. Finally, 
there is another aspect which is relevant in this context. When non-perturbative field 
configurations (domain walls, skyrmions, sphalerons, etc.) play a role, these can be taken 
into account with inhomogeneous background fields (i.e. mean field realizations). This 
may also be important for thermalization. 

In our previous work [17] the initial state was such that only a few of the low momenta 



modes of the mean field realizations carried all the energy. Such initial conditions were 
useful for equilibration tests. We found partial thermalization to an approximate Bose- 
Einstein (BE) distribution, for a limited time. Gradually the distribution also started to 
show classical equipartition features. For practical applications it is therefore important 
to get more information on the relevant time scales, and this is the main subject of this 
paper. 

In the present work the initial energy is distributed over all momentum modes of the 
realizations as in (H). This initial density operator is still out of equilibrium because of 
interactions. As will be shown in the next sections, this leads to creation of quantum 
particles of all energies with a Bose-Einstein distribution. Depending on the interaction 
strength and temperature, this equilibrium may last for a long time before classical features 
begin to dominate, and we are able to better quantify the relevant time scales. 

An important diagnostic in this and our previous work is the particle distribution 
function, which is defined in terms of equal-time two-point functions. The definition 
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assumes a quasi-particle picture and we check this here by performing a Monte Carlo 
computation in imaginary time to compute the exact two-point correlation function. 

The computations of the quantum mode dynamics is numerically very expensive. 
Therefore we also study the effect of reducing the number of quantum field modes. With 
the present initial conditions in terms of a temperature T = it is natural to try 
eliminating modes with |k| T. The often used classical approximation corresponds to 
the extreme of leaving out all quantum modes and to see how close this can mimic the 
quantum world we also make a comparison with this case. It turns out that even with a 
substantially reduced number of modes, this extended Hartree method fares much better 
than the classical approximation. 

We again use the (p 4 model in 1 + 1 dimension as a test case. This model, the Hartree 
ensemble approach and initial conditions are briefly recalled in Sect. |2|. In Sect. we recall 
our definition of the particle distribution in terms of equal time correlation functions and 
perform a Monte Carlo check on the underlying quasi-particle picture. In Sects. |] and [5] 
we study equilibration time scales at weak and stronger coupling. Next, in Sect. ||, we 
show that the expensive computations of the quantum mode dynamics can be substan- 
tially reduced by using only a limited number of mode functions. In Sect. |7| we make a 
comparison with classical dynamics. Finally in Sect. |8] we discuss the results. 



2 Method 



2.1 Hartree ensemble approximation 

Consider the Hamiltonian of a scalar field in one dimension, discretized on a lattice, 

H = £[3*2 - + + JA#], (3) 

X 

with x = a, 2a, . . . , Na, A(p x = (ip x+a + x - a — 2ip x )/a 2 . The volume is L = Na with 
periodic boundary conditions; The Heisenberg equations follow as, 

fix = Kx, fix = (A<£) x - fi 2 x ~ >4>x- ( 4 ) 

Rather than solving the exact operators from the Heisenberg equation, we use the gaus- 



sian or Hartree approximation (see e.g. Wll for details). This amounts to approximating 



the field operators as linear combinations of time- independent creation and annihilation 
operators 6* and b a , 

<px = ^+£[/^+/r%L 

a 

fix = ^+Yl{ffba+frbi], (5) 

a 

with time-dependent mean fields </? and tt, and mode functions f a . All information is 
contained in the one- and two-point functions, which can be expressed as 

(0x) = <p a , C xy = (0 x y ) - (0 x )(0 y ) = £[(1 + nDfiff + nlftf?], (6) 
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and similarly for the one- and two-point functions involving n. The n° a are the particle 
number densities in the initial state. In our numerical work we shall take n° = 0, such that 
the system is described by a pure state wave- functional. All 1PI higher-point functions 
are zero (i.e. higher-point functions factorize in one- and two-point functions). The mode 
functions represent the width of the wave-functional, allowing for quantum fluctuations 
around the mean field. Alternatively one can think of them as describing the particles 
in the model. Since the Hartree method uses (gaussian) wave-functionals, we expect to 
improve on the classical dynamics. Of course we cannot expect to capture all quantum 
effects, e.g. tunneling is beyond the scope of this gaussian approximation. 

Substituting the gaussian ansatz (Q) in the Heisenberg equations @ and taking ex- 
pectation values, we find self-consistent equations for the mean field <p x and the mode 
functions 



x 



A Vx - /i 2 + A^, + 3A]T(2< + 1)|/«| 2 



X ■ 



't = A/-- [// + 3A^ + 3A^(2n° + l)|/:| 2 ]/« (7) 

a 

To solve these equations we use a leap-frog algorithm, which is stable for sufficiently small 
time steps. Since there are N mode functions (the lattice has N sites), the amount of 
work to progress the fields over one time-step is 0(N 2 ). 

Above we described the Hartree approximation. In the Hartree ensemble method this 
approximation is applied to each individual realization \ip, ir)((p, 7r| of the initial conditions 
as in ([!]). So in eq. (ffl) the gaussian brackets stand for (•) = (<p, ir\ ■ \<p, it) and the average 
over <f, 7r is only taken in the evaluation of observables. Furthermore these states are pure, 
hence the initial particle density n° = in (E|) and ([7|). 

In this way we compute correlation functions with a generally non-gaussian density 
operator p = X^Pil^ j 71 " )(v >7T^|, as 

Sxy = $>[cS + <P { Sf\ - (E^) (E^'O- (8) 

i i j 

The C X y and ip x 1 ^ are computed with gaussian pure states using @. This means that 
in the time-evolution the gaussian approximation is used, while expectation values are 
calculated using the more general initial density operator. 

It should be stressed that for typical realizations the mean field (fx is inhomogeneous 
in space, in contrast to the ensemble average £\ Pi<p x ^ which is in fact homogeneous for the 
initial conditions we shall employ. We also note that the equations (0) can be derived from 
a hamiltonian. Since the equations are also strongly non-linear, this might suggest that 
the system could evolve to an equilibrium distribution with equipartition of energy, as in 
classical statistical physics. On the other hand, the mode functions satisfy Klein-Gordon 
type orthogonality and completeness relations that could obstruct such an equipartition 



and it is not easy to predict the equilibrium distribution of the model |1 1 
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2.2 Initial conditions 

In order to solve the equations of motion ([/]), we must specify initial conditions for the 
mean field <pW and the modes f a ^ of the individual Hartree trajectories as well as the 
weights Pi. This amounts to specifying the initial density operator p. As explained in the 
introduction, we use coherent (pure) states to represent the initial density operator, hence 
we use initial modes functions as in the vacuum state (with the initial particle density n° 
equal to zero). The initial mode functions may be taken as plane waves with wave vector 
k (a — > k), with a normalization appropriate for a zero temperature free field, 

lt(o) = f k M = -*"*-7s=r (9) 



(the same for each realization i). L is the system size and tu k = ^Jm 2 + [2 — 2 cos(aA;)]/a 2 , 
with m the zero temperature mass. We will not choose the initial p as far out of equilibrium 



as we did in ref. | 17| . There we took a superposition of only a few low momenta modes, 



J max 



= AmJ2 cos(2ttjx/L - ifjf), (10) 



<py = V 



where v = ^m 2 /2A is the vacuum expectation value of the mean field in the "broken 

symmetry phase" at zero temperature and ip^ is a random phase. Here we choose the 
mean fields from an ensemble with a Bose- Einstein distribution for the tp and 7r momentum 
modes as in @, 

p k {p k , 7T k ) oc exp[-(e^ /To - 1)K 2 + uft<p k - 5 k0 v) 2 )/2u k }. (11) 

Then the initial density operator is that of a free field thermal quantum ensemble (cf. 
appendix A in fl7H), 

P = II / d Vkdir k p k {p k , 7T k )\(p k , 7r k )(p k , n k \ oc e~ H ° /T ° . (12) 
k 

It should be stressed that this ensemble is not in equilibrium, even though the particle 
densities we compute from the initial conditions (after averaging over a large number 
of realizations) have a BE distribution. This is clear, because the mode functions do 
not contribute at all to the initial particle density. In each individual run we therefore 
expect quick excitation of the mode functions from their vacuum state, i.e. quantum 
particles will be created, using energy from the mean field. Moreover we use the free 
field dispersion relation uj\ = m 2 + [2 — 2 cos(ak)]/a 2 in the initial distribution (|TT|), with 
the zero temperature mass m = m(0). In thermal equilibrium this should become the 
temperature dependent mass m(T) of the quasi-particles. Nonetheless we expect that 
these initial conditions will lead to a much faster thermalization than initial conditions of 
the form (|T0|). 
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Figure 1: Dispersion relation computed from a Monte Carlo simulation of the euclidean 
time version of the model. The model parameters are: X/m 2 = l/2v 2 = 1/4, Lm = 25.6, 
1/am = 10 and T/m = 1, with 20 steps in the euclidean time direction. Here k is 
the lattice momentum \/2 — 2 cos(ak)/a. The statistical error bars are smaller than the 
symbols. The right figure is a zoom-in at small k where deviations from linear behavior 
are visible. 



3 Monte Carlo check 

To study time-dependence of particle energies and densities we define energies u k and 
densities n k from the Fourier components of the (p and % two-point functions, 

jJ2 e ~ lk(x ~ v)s -y = K + i)M- 

xy 

}_Y^ e -iH*-y) Uxy = [n k + \)u k , (13) 

xy 

where U xy is similarly defined as S xy in (||), with replaced by ir. For weakly coupled fields 
one expects the equilibrium particle densities defined in this way to have a Bose-Einstein 
distribution, and energies to have a free quasi-particle dispersion relation, approximately, 

n k = ( e ^/T _ ^2 = m(T)2 + [ 2 _ 2 cos(ak)}/a 2 . (14) 

The effective mass m(T) of the quasi-particles is temperature dependent. In the following 
we shall use the zero temperature mass m = m(T = 0) to scale dimensionful quantities. 

To substantiate this expectation (|T4"D, we have performed several Monte Carlo simula- 
tions of the euclidean time version of our model at parameter values in the same range as 
we use for the Hartree simulations. In Fig. [I] we show the dispersion relation computed 
from such a Monte Carlo simulation. We chose a temperature T/m = 1 and measured 
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Figure 2: Particle densities as a function of energy, plotted as log(l + 1/n). In Figs, a-c 
on the left the initial T /m = 1; on the right the initial temperature is high, T /m = 5. 
The model parameters are: X/m 2 = l/2v 2 = 1/12, Lm = 25.6, 1/am = 10. 



S xy . We stress that such a Monte Carlo simulation gives the exact (up to statistical errors) 
results for the finite temperature Green function. Making the assumption that rik has the 
BE form, we computed the ujk from S xy using (0). As can be seen from the figure, the free 
form ([14]) for the quasi-particle dispersion relation holds very well, with m(T)/m m 0.43. 



This value is close to that found with the effective potential calculations in the Hartree 
approximation [17J, which gives m(T)/m ~ 0.41 at T/m = 1 (see also Fig. ||). The effects 
of the temperature and interactions show up almost exclusively in the value of the effective 
mass m(T). 
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4 Weak coupling 



In our previous work, using the far out-of-equilibrium initial conditions (0), we found 
that particles of increasingly higher energy are created and acquire densities with a BE 
distribution. However, this thermalization progressed rather slowly to high energies, such 
that the low momentum particle densities already started to deviate from a BE distribution 
before particles with energies of a few times the temperature could participate in the 
equilibrium. These two phenomena - particles being created with densities that have a BE 
distribution and the gradual emerging of equipartition-like features - will be investigated 
below using the thermal initial conditions (1121) . 

To probe the large time behavior we shall use stronger coupling and higher energy 
densities than in our previous work. However, first we show results at the same coupling 
as used in |Tj]]. The coupling constant X/m 2 = \/2v 2 = 1/12 is in the "broken symmetry 
phase" of the model. The volume Lm = 25.6 and the lattice cut-off 1/am = 10. 

We plot log(l + l/n.fc) rather than rik itself, because in this way a BE distribution shows 
up as a straight line with a slope equal to the inverse temperature. The scattering in the 
data points is due to using only a few Hartree realizations (only two initial conditions). 
At low temperature, T /m = 1, Figs. |2]a-c, the evolution is very slow and there is hardly 
a sign of emerging classical features even at the largest time tm pa 50000. Even though 
the particle distribution does not change, there is a persistent, slow transfer of energy 
from the mean field into the modes. At tm = 200, 50% of the energy is still in the mean 
field, at tm = 6000 this has dropped to 25% and at tm = 50000 it is still some 15%. 
The effective mass stays roughly constant, m(T)/m pa 0.94, which is consistent with the 
effective potential for To/m = 1. 

At higher temperature, To/m = 5, but with the same weak coupling, there is again a 
wide window in which the particles have a BE distribution without significant distortions 
(Figs. 0d-f). However, we see classical-like features emerging for tm > 4000: compared 
to the BE distribution, the low momenta modes become under-occupied, while the high 
momenta modes become over-occupied. We find that at the latest time tm pa 50000 the 
distribution for u/m < 7 can be described reasonably well with an ansatz = co+Td/cui-. 
Without the constant cq pa 0.25 the fit would be poor. 

In this simulation we find an interesting behavior of the effective mass, shown in Fig. |4j. 
For comparison, we also show in Fig. [3] the effective mass calculated using the Hartree 
effective potential at the same model parameters. First the mass is steadily decreasing, 
which is appropriate when the temperature is decreasing and the system is in the hot, 
symmetric phase. At tm pa 14000 there is a sharp turnover and the mass starts to increase 
as in the cold, broken phase. The temperature at that point T c i/m pa 1.6, obtained from 
a classical fit, is close to the temperature T c /m = 1.8 of the first order phase transition 
computed from the effective potential. [[] Also the average mean field fluctuates around zero 
before and around v pa 1.8 after the transition, reasonably close to the effective potential 
prediction v pa 2 for T/m ~ 1.6 — 1.8. The reasonable quantitative agreement between the 
simulation, which shows classical features, and the effective potential computation, which 

■ We recall that in the exact theory there would be a cross-over instead of a first order phase transition. 
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Figure 3: Temperature dependence of the effective mass computed using the Hartree ef- 
fective potential fcF^j, X/m 2 = 1/4 (solid line) and 1/12 (dotted line), mL = 16 (volume 
dependence is very small). 
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Figure 4: Time dependence of the effective mass m(T) for the same model as shown 
in Figs. ^d-f. The mass is determined as the lowest energy ujq (dotted line) or from a 
quadratic fit to the dispersion relation (full line). 



assumes a BE distribution, illustrates that the thermal mass is dominated by the low- 
energy particles, for which there is little difference between a BE and classical distribution. 

5 Stronger coupling 

We now turn to the stronger coupling X/m 2 = l/2v 2 = 1/4, in order to make processes 
evolve faster. In Fig. |5|a we show particle densities nk computed only from the mode 
functions. We ignore the contribution from the mean field in (|8|), because we want to 
focus on the particles described by the mode functions. In Fig. ||a one sees that already 
after a short time, tm > 10, particles have been created over a wide range of energies, 
uj/m < 6. The densities are reasonably well described by a BE distribution with a time 
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Figure 5: The same as Fig. [|, but at a stronger coupling X/m 2 = 1/4. 

dependent temperature. This temperature initially increases rapidly from T/m = at 
tm = 0toT/m^0.6at tm = 10 and then gradually increases further to T/m pa 0.9 at 
tm = 100. (Recall that the temperature of the initial condition is Tq = m.) 

Figs. P>c, which are obtained using both the modes and mean fields to compute the 
correlation functions, show that the densities of particles with large momenta tend to 
remain at a BE distribution also for later times, with a very slowly increasing temperature 
T/m = 0.93 — 1.13. However, one also clearly sees deviations from the BE distribution 
developing, starting at the low a;-side of the spectrum. 

From these data we infer two time scales: First there is the rate at which the tem- 
perature of the BE distribution of the quantum particles is established. Second there is a 
rate at which the classical-like distribution sets in. Fig. ^| shows the time dependence of 
these two processes. The BE temperature was computed by fitting log(l + 1/n) = uj/T 
(only using the mode function contribution) for 2 < u/m < 4. The classical temper- 
ature was found from fitting n = T c i/uj for u/m < 2. The time dependence of these 
temperatures is reasonably well described by an exponential approach to an equilibrium 
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Figure 6: Time dependence of BE (from the modes only) and classical (from the modes 
and mean field) temperatures for the data of Figs. \^a-c. 



value, T BE (t) = A - Be~ t / TBE and T d (t) = A' + B'e'^. We find mr B E ~ 20 and 
mT c i pa 2500, showing quantitatively that the approximate BE thermalization happens 
much faster than the emergence of classical-like behavior (Note that T c i becomes much 
lower than Tbe, which itself is somewhat smaller than To, in agreement with the eventually 
expected classical equipartition). 

At higher temperature, the distribution roughly follows the same development. Sur- 
prisingly enough the distribution keeps its approximate BE form much longer, while at a 
higher temperature one expects a stronger effective coupling, and thus shorter time-scales. 
In Figs. [|d-f the initial temperature is T Q /m = 5. At this higher temperature and on 
a correspondingly larger energy scale, the deviations from a BE distribution appear less 
pronounced at early times. But even at tm = 4600 the particle densities are reasonably 
well described by a BE distribution with a temperature T/m pa 4.8. At this time, there is 
a small reduction of the density of low momentum particles (n is up to 15% smaller than 
the BE density, which is hard to see on the log-plot). At the same time the density of 
particles, with uo/m in the region 10 — 12, increases a little. This trend continues and at 
tm = 22600 there is classical behavior for uj/m < 12. 

The dashed line in Fig. [5|f is a fit of the form n = T c i/u, which gives a "classical" 
temperature T c i/m pa 2.1. The good quality of this fit for u/m < 12 suggests that the 
BE distribution gradually turns over into classical equipartition. However, for still larger 
times the distribution is no longer well described by a simple n oc 1 / cj-dependence. We did 
not determine the final equilibrium distribution, because of the extremely long (computer) 
time this would require. 

In Table 1 we summarize our results for tbe and t c i, including also fits to the time 
dependence of the particle density Ylk n k> computed from the modes only or mean field 
plus modes, as in Figs. || These results do not show a clear dependence on the coupling or 
temperature, contrary to the expectation of much smaller time scales at higher temperature 
and/or stronger coupling. We believe that this is accidental, due to the fact that at 
T /m = 5 and/or X/m 2 = 1/4 the system is in the "symmetric phase", while it is in 
the "broken phase" for T /m = 1 and X/m 2 = 1/12. We have noticed previously [fT7|| 
that in the "symmetric phase" the system evolves much more slowly than in the "broken 
phase". Note that the numbers in the table are subject to systematic uncertainty, since 
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X/m 2 


1/12 


1/4 


To/m 


1 


5 


1 


5 


rriT BE 


35 


35 


25 


25 


mTd 


> 15000 


3000-5000 


2500-3500 


2000-5000 



Table 1: Results for the BE equilibration time, Tbe, an d the time scale for the drift towards 
classical equipartition, r c i, obtained from fits to Tbe and T c i as in Fig. as well as similar 
fits to the time dependence of ^ fc nk ■ 



the mode system starts far from equilibrium and the time dependence not always follows 
an unambiguous exponential relaxation. This applies in particular for the simulation at 
X/m 2 = 1/4, To/m = l, which is very close to the "phase transition". 

Besides looking at the particle number distribution, it is interesting to follow the ef- 
fective mass m(T) in time. Comparing it with the temperature dependence computed 
analytically using the Hartree effective potential fIT|| , gives another measure for the ef- 
fective temperature of the system. The simulation of Figs. |5]a-c gave an effective mass 
which increased slightly in the range m(T)/m = 0.84 — 0.89. From the effective potential 
calculation we then infer that the temperature should be in the range 0.5 < T/m < 0.7, 
i.e. in the "low temperature" phase of the model, cf. Fig. [3], which is confirmed by check- 
ing the values of the mean field. This temperature is considerably lower than the BE 
temperature T/m = 1.0(1) estimated from the particle distribution at higher momenta, 
but is consistent with the temperature obtained from fitting rik at the smaller Uk with a 
classical distribution. The same is found for the high-temperature simulation of Figs. ||d- 
f: m(T)/m decreases from 1.12 at the start to 0.60 at tm = 22600, which corresponds, 
using the effective potential, to a decrease from T/m m 5 to T/m ~ 2, consistent with 
the observed T c i/m m 5 to T c [/m « 2.1. As mentioned before, the difference between the 
classical and BE distribution is unimportant for the dominant nk, those at low momenta. 

6 Reduced number of mode functions 

If the positive results for the performance of the Hartree ensemble method at shorter times 
carry over to more realistic models in 3 + 1 dimensions, one has to confront the problem 
of the high computational cost of this approach. Taking the continuum limit on a finite 
volume in d dimensions, i.e. increasing the number of lattice sites N in each direction, the 
cost of our approach increases oc N 2d+1 : There are 0(N d ) fields which have to be updated 
O(N) times, assuming a fixed value of the time-step ao/a. 

Most of this cost comes from having to solve all N d mode functions f a . However, 
many of these modes would represent particles with very high momenta \k\ ^> T. Such 
particles have very low densities and should be irrelevant for the physics at lower scales. 
This suggests reducing the number of mode functions in our simulations. We have tested 
this idea by comparing a simulation on a lattice with iV = 128 sites using all 128 mode 
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Figure 7: Particle densities at tm = 50, 90, 300, obtained from simulations using the full 
number of modes (drawn lines) and using only modes for which lo^/m < 17 ~ 3T/m. The 
left figure shows log(l + 1/n), the right shows the density n itself (Lm = 5.7, 1/am = 22.3 
and X/m 2 = 1/12). 

functions, with the same simulation using only 32 mode functions. This induces a max- 
imum energy uj mSLX /m pa ^2 — 2 cos(327r/128)/am pa 17, which is much larger than the 
temperature T/m pa 6.8 that we will use in the test. The uj maK here refers to the energy 
of the initial mode functions which are plane waves. 

In order to make as detailed a comparison as possible, we show the results for the 
particle density obtained from the mode functions only, leaving out the contribution of 
the mean field in (H). As shown in Fig. [7], this partial distribution changes with time 
during thermalization (cf. Fig. |5]a). The drawn lines represent the data obtained with the 
full number of mode functions. The dots represent the results obtained using the reduced 
number of modes. The left figure shows the familiar log(l + 1/n) form of the density, the 
other figure shows the density n itself. As is evident, these results reproduce the data 
from the reference simulation accurately up to uo/m pa 12, which is close to uJ max /m pa 17. 
Notice that the densities, computed from eq. (pf), drop ton = —0.5 for u/m > 17 (Fig. [7] 
right). For these high momenta there are no more mode functions available to provide the 
vacuum fluctuations that should lift the density to zero. It indicates that at high momenta 
there is still a roughly one-to-one correspondence between mode functions and momentum 
labels of the particles. 

7 Classical approximation 

Even using fewer mode functions, the Hartree approach is much more expensive than the 
classical approximation (which has no mode functions). So it is important to check if our 
results cannot in some way be mimicked by a classical approximation. The standard way 
to implement the latter at high temperature, is to average over initial configurations drawn 
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Figure 8: Same as Figs, ^d-f, but using classical dynamics. Figs, a-c have BE-type initial 
conditions with temperature T /m = 5, in Figs, d-f only a few of the low momentum 
modes of the mean field carried all energy. The other model parameters are: X/m 2 = 1/4, 
Lm = 25.6 and 1/am = 10. 

from the Boltzmann distribution. Up to modifications by the interactions this implies a 
classical distribution function n{u) = T/u, with a slow fall off causing Rayleigh- Jeans type 
divergences. Actually, in 1 + 1 dimensions such divergences are absent in (/^-correlation 
functions 0. 

Here we want to ask a somewhat different question: to what extend can classical 
dynamics be used to represent a thermalized system with a Bose-Einstein distribution for 
the particle densities? To investigate this we shall use the same BE-type initial conditions 
(O) as in the Hartree case, as well as the much more out-of-equilibrium conditions of 
the form ( [10]) . We perform similar simulations and analyses as before, but now without 
mode functions (and with + 1/2 — >• nk in (|l~3|), as there is now no quantum vacuum 
contribution). Typically, the classical dynamics produces data with more noise, since the 
average contribution of the many mode functions tends to smoothen results in the Hartree- 
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dynamics. We counter this noise by averaging over 40-50 initial conditions, which is much 
more than we typically use with Hartree dynamics. We note in passing that the necessity 
to use a larger ensemble to obtain data of the same quality as with Hartree dynamics, 
diminishes the computational advantage of using classical dynamics considerably. 

At the same weak coupling and low temperature as in Figs. |2]a-c, we find that the 
classical system also preserves the Bose-Einstein distribution of the initial conditions very 
well. Even at the largest simulated time, tm = 50000, there is no compelling sign of 
equipartition in the particle distribution. This, however, may only show that the relaxation 
time-scale of the classical dynamics is very long, cf. ||, and that we are seeing remnants 
of the initial condition rather than thermalization. 

To speed-up the dynamics we increased the temperature to To/m = 5 and the coupling 
to X/m 2 = 1/4, as was used in Figs. |5|d-f. The results are shown in Figs. §a-c. Now 
the initial BE-distribution still persists for some time, but already at tm « 600 there are 
clear signs of equipartition setting in, whereas effects of similar magnitude only emerged 
at tm > 6000 with Hartree dynamics. The gradual move from the initial state towards 
classical equipartition happens much faster than in the Hartree ensemble simulations. 

Of course one might argue that this initial persisting of the BE distribution is of little 
significance, since it only demonstrates that it takes time to loose the effect of the initial 
conditions. In more realistic models we might not be able to specify initial conditions 
sufficiently close to thermal equilibrium and then one may not expect to encounter a BE 
distribution. Yet, somewhat surprisingly, starting with the far out-of-equilibrium initial 
condition (ITUI), we see that as the model steadily moves towards classical equipartition, 
the particles are distributed in a BE-like way in an intermediate stage. This can be seen 
in Figs. |8]d-f, where we show particle density distributions at simulation times in the range 
tm = 200 — 20000. At tm around 800 the particle densities follow a BE-like distribution 
over a wide range of energies. The coupling strength in this simulation is X/m 2 = 1/4, at 
weaker coupling this intermediate stage with BE-like distribution persists longer. However, 
it always smoothly turns towards classical equipartition on much shorter time-scales than 
when using Hartree dynamics (although, as mentioned above, the final equilibration to 
the classical distribution takes place on a very long time scale). 



8 Discussion 

Using Monte Carlo methods we checked that the quasi-particle assumption, which is basic 
to our definition of the particle distribution function, is well satisfied in the range of 



couplings and temperatures considered here and in |17|. 

From the results in this work combined with |T7| , the following picture has emerged for 
the 1+1 dimensional Xip A model in the "broken phase". The initial energy, which is put 
solely in the mean field of a realization, is subsequently transfered to the mode functions. 
This process takes place fairly locally in momentum space, i.e. mean field modes with 
momentum k excite primarily particle modes with momenta close to k, and the modes 
then thermalize locally to a BE distribution. In our previous work this approximate 
thermalization was more conspicuous because the initial distribution was further out of 
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equilibrium. Here the BE distribution was put into the initial condition for the mean 
fields. However, the corresponding density operator is still out of equilibrium because of 
the "wrong" initial thermal mass. The thermalization process is fairly rapid, within a 
time tbe ~ 25 — 35m _1 , for A/m 2 = 1/4, 1/12 and Tq/tti = 5, 1, as determined from the 
time-dependence of the BE temperature or the particle density J2k Uk ( c ^- Table 1). 

This time scale is similar to our findings with initial mean fields containing only low 
momenta [ I7[ . The subsequent thermalization of higher momenta is very slow. We ascribe 
this to a weakening of the non-linearities when the mean field looses much of its energy. 
When the mean field fluctuates around its (temperature dependent) equilibrium value, 
with diminishing amplitude, the dynamics becomes approximately that of Hartree with 
a homogeneous mean field, suggesting lack of thermalization. This also explains why 
the evolution to a classical-like distribution is much slower with the Hartree ensemble 
approximation than using classical dynamics. 

However, the fluctuations die out very slowly and even at very large times of order 
10 4 m _1 there is still 0(10%) of the energy in the fluctuating mean field. Nonlinear fluc- 
tuations remain, which lead eventually to classical-like equipartition (according to the 



effective hamiltonian and conserved "charges" |[17 



The time scale for such classical equipartition setting in could not be determined in 
17] , its determination is one of the results of the present work. We find that the system 
remains in an approximate quantal thermal state for times of the order r c / > 100 tbe (cf. 
Table 1). 

This is an encouraging result. For example, in a crude application of our 1+1 dimen- 
sional results to 3+1 dimensional heavy ion collisions, identifying m with the mass of the 
a-resonance m a = 600 - 1200 MeV, say 900 MeV, a time-span of 100t be = 2000m" 1 
would correspond to a reasonable length of about 450 Fermi. Within such a time-span 
the Hartree ensemble method may be a definite improvement on the classical dynamics 
usually employed for e.g. the "disoriented chiral condensate" . 

For application to 3+1 dimensions it is important that the numerical efficiency of the 
Hartree method can be significantly improved by using only a limited number of mode 
functions, corresponding to particles with sufficiently high densities (see sect. |6|). 

Leaving out the mode functions altogether, i.e. using classical dynamics, the results 
were qualitatively similar to those with Hartree dynamics, but the emergence of classical 
particle distributions goes faster by roughly an order of magnitude. So this may not be 
good enough for practical applications. 

With respect to thermalization it is good to keep in mind that in the Boltzmann 
approximation, the collision term corresponding to 2 — 2 scattering is identically zero, 
due to kinematical constraints in the ip 41 model in 1+1 dimensions. So thermalization has 
to come from inelastic scattering and/or off-shell effects. It is then important to realize 
that such effects are more pronounced in the "broken phase" of the model, which has a 
three point vertex and finite (as opposed to zero) range interactions. As mentioned in 
||17|| , thermalization is drastically less efficient in the "symmetric phase" at similar values 
of A/m 2 . It is sobering to recall the huge thermalization times found in || in the classical 
approximation, in the "symmetric phase". For example, using A/m 2 = 1/4, T/m = 0.2, 



16 



the formula (rewritten in our conventions) l/mr c i ass = 5.8 10 _6 (6AT/m 3 ) 1 ' 39 found in this 
work leads to a relaxation time mr c i ass « 1.3 10 5 . This is much larger than the mr c / pa 2500 
found here in the "broken phase" (Sect. 0). 

An interesting question is how the two time scales tbe and t c \ are related to particle 
scattering and damping. A perturbative computation (which includes direct scattering 
through the setting-sun diagram), indicates that the damping time would be of the order 
of the BE-relaxation time tbe (i-e. much shorter than the relaxation time away from BE 
behavior). Preliminary numerical results for the damping time are consistent with these 
values 0,0. 

This is a favorable result for the Hartree ensemble method. However the gradual 
drift away from a BE distribution and the corresponding cooling of the system reveals a 
shortcoming. This is additional to the incorrect prediction by the Hartree method, of the 
order of phase transitions. Further improvements are needed, in particular if large time 
scales are to be investigated. 
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